Experiments with a Malkus-Lorenz water wheel: Chaos and Synchronization 
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We describe a simple experimental implementation of the Malkus-Lorenz water wheel. We demon- 
strate that both chaotic and periodic behavior is found as wheel parameters are changed in agreement 
with predictions from the Lorenz model. We furthermore show that when the measured angular 
velocity of our water wheel is used as an input signal to a computer model implementing the Lorenz 
equations, high quality chaos synchronization of the model and the water wheel is achieved. This 
indicates that the Lorenz equations provide a good description of the water wheel dynamics. 
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I. INTRODUCTION 

Since the discovery by Lorenz"'^ that a simple three- 
variable set of ordinary differential equations can give 
rise to exceedingly complex behavior, the study of chaos 
has continued unabated. Over the years, chaos has been 
found in a variety of naturally occurring systems and 
many fascinating table top experiments have been con- 
ducted in order to quantitatively study aspects of chaos. 
Experiments range from dripping faucets^ and pendula"^ 
to chemical reactions'' and lasers.^ 

For purposes of practical applications of chaos, opti- 
cal and electronic systems operating on time-scales of 
nanoseconds or less are a current research focus. Yet 
for purposes of gaining intuition, chaotic mechanical sys- 
tems operating on timescales of seconds are unsurpassed 
because of their palpable mechanisms and the direct ex- 
perience they provide. A prime example of such an ex- 
periment is the Malkus-Lorenz water wheel, slowly rotat- 
ing one way and then another in an at once calming yet 
interestingly unpredictable fashion. 

This water wheel was first envisioned and constructed 
by Malkus and coworkers as a mechanical analogue of 
the Lorenz equations. ^^^^ It has fascinated many students 
(and teachers) and has become particularly well known 
since its discussion in Steven Strogatz's introductory text 
on nonlinear dynamics. The Malkus water wheel de- 
scribed in this paper also grew from such inspiration, 
starting as a senior thesis project of one of the authors 
(Rachel Fordyce). 

The most basic Malkus water wheels consist of a few 
leaking cups attached to the rim of a freely turning wheel 
whose axis may be either horizontal or tilted. A single 
stream of water at the top of the wheel will add water 
whenever a cup is close to the top. Although the mo- 
tion of such simple water wheels can be analyzed, their 
dynamics differs from that of the Lorenz model, and we 
refer to them as non-ideal. The goal in our experiment 
was to construct an ideal water wheel, one whose dynam- 
ics is described by the Lorenz equations. 

In this paper we aim to evaluate how close we have 
come to this goal. We present in Sec. II details about 
our experiment, with a focus on confirming three crucial 
assumptions that are made in deriving the Lorenz model 




FIG. 1: (Color online) Experimental implementation of 
the Malkus-Lorenz water wheel, constructed using a bicycle 
wheel, syringes, and rare earth magnets forming a magnetic 
brake. The schematics on the right show a side view (bottom) 
and a top view (top) of the wheel. (See text for details.) 



for the wheel's dynamics. The derivation of the model 
and a brief overview of some of its properties are pre- 
sented in Sec. III. Experiments resulting in both chaos 
and periodic oscillations are shown in Sec. IV. In an effort 
to make a more precise statement regarding the match of 
model and experiment, we turn in Sec. V to an analytic 
proof and experimental demonstration of chaos synchro- 
nization. We conclude with a discussion in Sec. VI. 



II. EXPERIMENTAL SETUP 

In designing our water wheel we strove to balance the 
desire for a simple and inexpensive experiment and our 
goal to come as close as possible to an ideal water wheel. 
Three properties of an ideal water wheel require consid- 
eration. First, the amount of water entering the cups 
per unit time has to be constant, which means that gaps 
between the cups and cup overflow have to be avoided. 
Second, damping is assumed to be entirely due to fric- 
tional torque proportional to velocity (viscous damping) , 
which means that dry friction (kinetic and static friction) 
has to be negligible in the experiment. And third, the 
cup leakages should be proportional to the water mass in 
a cup, which suggests that the outflowing water should 
exit through a pipe with a laminar flow described by the 
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Hagen-Poiseuille equation. With these requirements in 
mind, a water wheel was constructed as described below. 

As a wheel, a smallish 16.5 inch BMX bicycle wheel was 
used. It was chosen for its smooth bearings that guaran- 
teed very low axle friction. The wheel was mounted to a 
pair of hinged boards, enabling us to vary the inclination 
angle a of the wheel with respect to the horizontal (see 
Fig. 1). 

Fifty-six syringes of 50 cm^ volume each served as leak- 
ing "cups" and were attached to the perimeter of the 
wheel using a metal synching band. For each syringe, 
the sides of the finger pull was shaved off in order to al- 
low a snug fit with no gap in between the syringes. Sub- 
sequently to attaching the syringes, the wheel was bal- 
anced with small weights, such that the angle at which an 
empty tilted wheel would come to a stop was distributed 
uniformly. We found that rebalancing was necessary in 
order to remove the wheel's tendency to preferentially 
stop with its welding seam at the bottom. 

The motion of the wheel was detected by a rotary en- 
coder attached to the wheel's hub. Its output pulses were 
counted continuously by an interface circuit containing a 
HCTL-2016 quadrature decoder and counter chip that 
includes a 16-bit memory for storage. The counts on the 
memory chip were then acquired by a computer, with a 
10 Hz sampling rate for the data shown in this paper. 

To achieve viscous damping, a (nonmagnetic) Alu- 
minum disc was attached below the wheel such that it 
co-rotated with the wheel. A second non-rotating disc 
with several rare earth magnets glued to its surface was 
mounted opposite the first disc. The result of this ar- 
rangement is a magnetic brake. When the wheel is in 
motion, the stationary magnets induce eddy currents in 
the spinning disc which in turn produce magnetic fields 
that oppose changes in magnetic fiux. A frictional type 
torque develops that is directly proportional to the angu- 
lar velocity. -'^•^"^^ The gap between the discs, and therefore 
the braking coefficient, is tuneable in our setup. 

To check whether damping in our wheel can be mod- 
eled as torque linear in the velocity, we measured the 
angle as a function of time for a slowing wheel without 
water and fitted the results. To obtain a model for the 
fit, we assume that the only relevant torques are from 
viscous friction due to the magnets and kinetic friction. 
The magnetic brake contributes a torque Tvf — — kw, 
where k is the viscous friction parameter and a; is the 
angular velocity. "^^ The magnitude of the kinetic friction 
force is -Fkf — M , where /i is the friction coefficient and 
-F/v is the normal force. Since Fjq can be taken as con- 
stant and the setup is rotationally symmetric, the torque 
due to kinetic friction is a constant (|Tkt| — T'kf — const.) 
opposing the motion of the wheel. The evolution of the 
magnitude of the angular velocity is therefore described 
by the first-order ordinary differential equation 

wh -TT = -Kt^ - Sgn w) Tkf, (1) 

at 

where /wh is the moment of inertia of the wheel with 
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FIG. 2; (Color online) Measurements of angle versus time for 
an empty wheel are compared to fits that take into account 
damping due to viscous friction (2 parameter fit) and a com- 
bination of viscous and kinetic friction (3 parameter fit), (a) 
Angle-data and fits, (b) Residuals of fits and angular velocity 
Lo as determined from the 2 parameter fit. 

the attached empty syringes. Equation (1) can be inte- 
grated, without loss of generality, under the assumption 
that a;(t = 0) = > 0, yielding 

..(*) = -!M + (^„ + I^)e-(-/^-)*, t<t.top. (2) 

The magnitude of the angular velocity decreases in time 
and, once it reaches zero at time t = tstop, the wheel stops 
because the torques due to viscous and kinetic friction 
are zero and static friction implies that some finite min- 
imal torque has to be applied to set the wheel in motion 
again. Integration of Eq. (2) yields the experimentally 
measurable angle as a function of time (t < tstop), 

^^ = -^^+^Lo + — ) (i-e-t'^/^-)*), (3) 

where we have taken 9{0) = 0. A convenient parameter- 
ization is achieved by introducing the viscous damping 
rate 7 = n/Iwh and parameter Q = r^f/K, yielding 

e = -nt + ^^^{i-e-^'). (4) 

In Fig. 2 we show two fits of the data to Eq. (4): a 
three-parameter fit corresponding to Eq. (4) and a two- 
parameter fit that corresponds to Eq. (4) with damping 
entirely due to viscous friction, that is with Tkt = = 0. 
It is seen in Fig. 2a that both fits provide an excellent 
description of the data and that they are essentially in- 
distinguishable to the eye. The effect of kinetic friction 
becomes only noticeable when the wheel is turning very 
slowly. This is seen from the residuals shown in Fig. 2b, 
which scatter around a value of zero for the first four sec- 
onds when the angular velocity is large, indicating that 
both models provide a good fit, but show a non-random 
deviation from zero at later times when the angular ve- 
locity is small. In particular, the residuals of the last 
two seconds before the wheel comes to a stop show that 
the model with just viscous friction (two-parameter fit) 
sightly overestimates the angle. Despite this, these re- 
sults clearly indicate that damping is dominated by the 
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FIG. 3: Measurements of the water volume in the syringe as 
a function of time and the two parameter fit of Eq. (9). The 
dotted line indicates the value of VJ, (see text), (a) Syringe 
without needle (VoB = 4 cm^). Fit-result: Vb = 29 cm^, 
k = 0.10 s"^ (b) Syringe with needles (Voff = 18 cm^) of 16 
gauge (Fit-result: Vo = 29 cm^, k = 0.017 s"^) and 18 gauge 
(Fit-result: Vo = 29 cm^, k = 0.006 s"^). (c) Geometry of 
the syringe, consisting of a cylindrical main body (radius rs, 
height h), a tapered section, and either a nozzle or a needle 
as exit "pipe" (radius Vp, height Zp). The total height of the 
tapered section and exit "pipe" is £. 

magnetic brake and that our water wheel comes close 
to the ideal one where damping is entirely due to viscous 
friction. To keep the dynamic model simple, we therefore 
neglect kinetic friction, from now on, and set r^f = 0. 

The moment of inertia of the wheel can be determined 
by performing a second measurement and fit for a wheel 
with a known set of additional weights of total mass m 
attached symmetrically to the wheel's rim.^^ Using the 
resulting fitted parameter 7 = K/{Iwh + mR^), the pre- 
viously determined 7, and 

Iwh = mR^ ^ , , (5) 

7-7 

we found that our wheel has a moment of inertia of /^h ~ 
0.11 kg m^. 

For an ideal water wheel it is important that the inflow 
is symmetric with respect to the centerline that divides 
the wheel into left and right. In order to achieve this sym- 
metry and, at the same time, minimize the probability 
of cup overflow while still being able to achieve sufficient 
infiow rates, we used a system of eight individually ad- 
justable spigots.'^® The spigots were placed symmetrically 
with respect to the top, spanning angles 9 e [— ^^o.+^o] 
with 9q 26°. A uniform distribution of mass fiux, de- 
scribed by 

fo ee[-TT,eo) 

Qi(^)^{^ Oe [-00,00], (6) 
[0 0ei9o,n]. 

was ensured by verifying experimentally that for a sta- 
tionary wheel all eight top cups fill at identical rates. The 
total mass flux Qtot was measured using a digital flow 
meter'^^ and was tunable with a range of 0.03 - 0.09 kg/s. 

Finally, the ideal wheel should have a leakage rate that 
is proportional to the water mass. Due to the fact that we 



use syringes as "cups" , water leaks out either through the 
small cylindrical nozzles at the end of the syringes, which 
suggests a theoretical description in terms of a laminar 
pipe flow of a viscous fluid through the nozzle, or alterna- 
tively, needles can be attached resulting in an even longer 
"pipe" . Assuming that the laminar pipe flow of the ex- 
iting water is the dominant effect, one expects the total 
volume V of water in the syringe (radius rg = 1.08 cm) 
to change according to the Hagen-Poiseuille equation^^ 

dt ?>vpw \ Zp J 

Here, Tp is the exit "pipe" radius, is the density, and 
1/ the kinematic viscosity of water. The change AP of 
the modified pressure P over the exit pipe length Zp is 
assumed to be entirely due to gravity, implying AP — 
g pw {h+£) (see Fig. 3c). The total volume V — nr'^h+Vh 
consists of the water in the cylindrical part of the syringe 
and an additional approximately 2 cm'^ volume Vb at the 
bottom. Taken together, this suggests that the volume 
in a syringe behaves as 

^^-k{V + V,f,) (8) 

where VoS — tnri — Vb is a constant offset term and 
the rate fc is a parameter that we flt (nominally k ~ 
r^g/Siyr'^ Zp). The fit model is obtained by integrating 
Eq. (8), yielding 

V(i) = Voe-'=* + Voff(e-'=*-l). (9) 

This model is valid for water levels within the cylindrical 
part of the syringe, i.e. for V{t) > Vb. The drainage 
of the small bottom volume cannot be measured with 
our setup but we is expect it to deviate somewhat from 
model (9) due to the tapering of the syringe and capillary 
effects. However, since the bottom volume is very small, 
we neglect any potential small corrections and assume 
that Eq. (9) describes the leakage of the entire syringe. 

To check the validity of Eq. (9) , we determined the to- 
tal water volume as a function of time by using motion 
tracking software to analyze digital films of the declin- 
ing water level in a vertical syringe with and without 
attached needles. The data and fit results are shown in 
Fig. 3. It is seen that the fit is quite accurate already 
for the syringe without needle, as shown in Fig. 3a, and 
becomes even better when a needle is attached, as seen in 
Fig. 3b. Note that by attaching needles and by changing 
the gauge of the needles the leakage rate k can be varied 
by orders of magnitude. No needles were attached to the 
wheel for all the data shown in this paper. 

Writing Eq. (8) in terms of the total mass of water 
in all 56 syringes, Mtot = 56pw V, and introducing an 
effective mass flux Q^ff via 

QesiO) = Qid) - 56/(2^) k pw Voff, (10) 
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where Q{0) is given by Eq. (6), one finds that the total 
mass in the wheel evolves according to 



dMtot 
dt 



(11) 



with — Qtot — 56 Vos- Equation (11) shows that 
the total mass of water is conserved in the asymptotic 
limit because Mtot = Q^/k for t — )• oo. 



III. THE MODEL 

To make this paper self-contained, we recall here the 
well known derivation of the equations describing the dy- 
namics of the water wheel and briefly discuss the param- 
eter dependence of the model solutions. 



A. Derivation 

In our derivation we follow Strogatz and model the 

water as being distributed as a continuous ring around 
the rim of the wheel. Matson^^ showed that the same 
result can bo obtained when considering an ideal water 
wheel with a discrete set of cups. 

In the continuum approximation the mass distribution 
Tn{9,t) aroimd the wheel's rim is defined such that the 
mass between the angles 9i and 02 (fixed in the lab frame) 
is 



M{t)= / m{e,t)de. 



The change in time of m{0,t), given by 



dt 



86 



(12) 



(13) 



has three contributions: the constant inflow Q{9), de- 
scribed by Eq. (6); the leakage, which is proportional 
to the water mass and which also contributes a con- 
stant offset to the inflow, resulting in Qeff as described 
by Eq. (10); and a third term that takes into account 
the wheel's rotation. The third term expresses the fact 
that the mass density at angle and time t + At, i.e. 
m{9,t + At), is identical to the mass density at time t 
at angle 9 — A9, i.e. m{9 — A9,t), for a wheel without 
water inflow and leakage, rotating with angular velocity 
u) such that A9 = oj{t) At. 

The change in time of the angular velocity is due to 
the applied total torque, which has three components. 
The magnetic brake contributes Tvf = —KLj{t), as dis- 
cussed in Sec. II. The infalling water contributes a torque 
Tspin up ~ —Qtoi ^(t) because it enters the wheel a dis- 
tance R from the center with zero angular velocity and 
is spun up to an angular velocity of co before leaking out. 
Gravity provides the third torque because each infinites- 
imal mass element m{9, t) d9 on a wheel that is tilted at 



an angle a contributes 



c^'^grav = Rq sin(a) sin(^) m{9,t) d9. 



(14) 



Taken together, one obtains 
duj{t) 



dt 



- (K + QtotR^) oj{t) 
+ Rgsm{a) f m{9,t)sm{9) d9, (15) 

J —IT 



where we have taken /tot to be a constant, which is valid 
after an initial transient, i.e. after the total water mass 
has come exponentially close to its final constant value 
[see Eq. (11)]. 

Since m{9, t) is periodic in 9, one can expand this func- 
tion into Fourier-modes, 



n=l 



m{9, t) = ^ a„ {t) sin(n ^) + -f- ^ 6„ (i) cos(n6'). 

(16) 



Under the assumption that Q{9) is truly symmetric with 
respect to the wheel's center line, the effective inflow Qefi 
can be written as 



Qeff(^) = I + E ^) - S'^'^wX^off, (17) 



56. 



n=0 



with coefficients 

Qiot . r n \ 

Qn = smc(n^o) 



27r 



(n = 0,l,2,...), 



(18) 



for the case that Q{9) is given by Eq. (6). Substituting 
these series into the partial differential equation (13) and 
integro-diffcrcntial equation (15) and equating the coeffi- 
cients of each harmonic separately, one obtains an infinite 
set of ordinary differential equations that describes the 
evolution of the angular velocity uj{t) and the Fourier- 
mode amplitudes o„(i) and b„{t). Amazingly, the evo- 
lution of the angular velocity and amplitudes ai{t) and 
bi (t) is entirely decoupled from the evolution of all other 
amplitudes. Therefore, the problem is reduced to a three- 
dimensional system of ordinary-differential equations 



du! 

lit 
da I 
~dt 
dbi 
~dt 



K + Qtol R- 



u;{t) 



TT R (J sin a 



hot 

kai{t)+Lj{t) bi{t) 

kbi{t)-uj{t) ai{t) 



ai{t) 



(19) 



It turns out that Eq. (19) can be mapped onto the Lorenz 
equations by introducing the dimcnsionlcss parameters 



Ik- 

k^ 



Qtot R? 

hot 



qi TT Rg sin a 
¥k + Qtot -R2 ' 



(20) 
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the dimensionless time s = kt, and the coordinates 

x=- j/= — fli z = p bi. (21) 

k qi qi 

In the new coordinates, Eq. (19) becomes 

X — a {y ~ x) 

y = px-y-xz, (22) 
z = X y — z 

where the overdot denotes the derivative with respect 
to the dimensionless time s. Equation (22) is identical 
to the Lorenz equations with the third parameter (often 
denoted as b) equal one. 

Above mapping of the problem onto Eq. (22) relies 
in an essential way on the fact that Eq. (13) is linear 
in m{0,t), explaining why an ideal water wheel has to 
have a leakage that is proportional to the water mass. 
Furthermore, to obtain the correct form of the Lorenz 
model, Eq. (15) needs to be linear with respect to uj{t), 
making a purely viscous damping a principal requirement 
of ideal water wheels. 



B. Characterization via Lyapunov Exponents 

Starting with Lorenz's seminal 1963 paper, ^ the Lorenz 
equations have been studied intensively over the last 
decades. What makes them so interesting is that this 
simple looking deterministic system has extremely rich 
dynamics. For some parameters the asymptotic solutions 
are steady state solutions, for other parameters periodic 
oscillations are found. -'^^'^^ In addition, and this is the 
reason that the Lorenz equations are famous, over wide 
parameter ranges the solutions exhibit a qualitatively dif- 
ferent and new form of behavior, they are chaotic. That 
is, they oscillate irregularly, never exactly repeating, and 
depend very sensitively on the choice of initial conditions 
(the initial values of the three variables) . Some rigorous 
results are available on the Lorenz equation's chaotic so- 
lutions and bifurcations, ^^"^^ where a bifurcation refers 
to a sudden qualitative change in behavior that occurs 
when a small smooth change is made to a parameter 
value. Yet, since such rigorous results are extremely chal- 
lenging to obtain, one generally resorts to approximate 
numerical calculations to characterize solutions. 

To numerically classify solutions of Eq. (22) as a func- 
tion of the parameters, we compute the largest Lyapunov 
exponents. Lyapunov exponents provide a useful mea- 
sure because they are independent of the particular initial 
conditions used in computing them and are invariant un- 
der smooth nonsingular coordinate transformations such 
as those in Eq. (21). In addition they directly measure 
a distinguishing characteristic of chaotic solutions, their 
"sensitivity to initial conditions". For chaotic systems, 
small differences in initial conditions yield after an expo- 
nentially short time widely diverging solutions, render- 
ing impossible the long-term prediction of outcomes from 
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FIG. 4: (Color) Shown is the numerically determined largest 
Lyapunov exponent as a function of p and a. White cor- 
responds to a zero exponent implying periodic oscillation of 
uj. Red colors correspond to positive and blue colors to neg- 
ative exponents, indicating, respectively, chaos and steady 
state evolution of ui. The approximate upper boundary of 
the experimentally accessible regions is shown by the solid 
black line. The yellow squares indicate parameter values cor- 
responding to the chaotic and periodic experimental time se- 
ries in Fig. 5. 

finite-precision measurements of initial states. Roughly 
speaking, two trajectories in phase space with initial sep- 
aration (5xo diverge as 

|<5x(OI«e^Vxo|, (23) 

and A, the largest Lyapunov exponent (LLE), is a posi- 
tive real number for a chaotic system. In contrast, sta- 
ble steady state solutions have a negative LLE, because 
perturbations away from the steady state decay expo- 
nentially, and stable periodic solutions have a LLE with 
value zero, because perturbations along a trajectory nei- 
ther grow nor shrink in time (and perturbations perpen- 
dicular to the trajectory decay exponentially). 

Figure 4 depicts the results of a numerical calculation 
of the LLEs, providing a map of the water wheel's dy- 
namic behavior as a function of the relevant effective pa- 
rameters, p and a. (Details about the numeric calcula- 
tion are given in the Appendix.) It is seen that for small 
values of either parameter only negative LLEs are found, 
corresponding to steady state behavior of the angular 
velocity. For larger parameter values, the wheel exhibits 
either chaotic or periodic behavior. We find that there 
exist large continuous regions with a zero LLE, shown in 
white, indicating periodic behavior. In addition, there 
are large continuous chaotic regions with positive LLEs, 
shown in red colors (dark gray). Yet, this numerical cal- 
culation also reveals that there are periodic windows em- 
bedded inside the chaotic region and that those windows 
seem to form a fractal set, with an increasing number 
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FIG. 5: (Color online) (a,b) The experimental time series of 
the rescaled angular velocity, x = uj/k, as a function of the 
dimensionless time s = kt is shown as the upper trace and 
the corresponding numerical solution Xnum as the lower trace. 
The parameters are (a) a ~ 3.6 and p ~ 140, resulting in 
periodic oscillations, and (b) cr ~ 2.5 and p ~ 66, resulting in 
chaotic oscillations. A corresponding two dimensional time- 
delay embedding of the experimental data is shown in (c) for 
the periodic and (d) for the chaotic case. 



of ever smaller periodic vifindows, each organized along a 
line in parameter space. This intricate structure means 
that the qualitative dynamics of the Lorenz system sen- 
sitively depends on the chosen parameter values. It also 
should be noted that for some parameter values several 
different asymptotic solutions (i.e. attractors) may co- 
exist in phase space. Our numeric result only depicts 
the LLE associated with one of those solutions, the one 
whose basin of attraction includes the randomly chosen 
initial condition of the simulation. 



IV. EXPERIMENTAL CHAOS 

Figure 4 is a useful guide for experiments because it 
informs our choice of wheel parameters that result in a 
desired dynamic behavior. This a-priori connection be- 
tween the wheel operating condition and expected dy- 
namic behavior is possible because we can explicitly write 
the effective parameters p and ct as a function of the 
experimentally tunable and fixed measurable parame- 
ters. Leaving aside the option of attaching needles to 
the wheel, the experimentally adjustable parameters are 



the inclination angle a, the viscous damping rate 7, and 
the total mass flux Qtot- We obtain p and a in terms 
of these quantities by substituting Eq. (18) for qi, writ- 
ing /tot as /tot = /wh + Mtot/?^ = /wh + Qcff ^Vfc with 
Qlg = Qtot - SefcpwV'oft, and using k = 7/wh- These 
substitutions yield 



p{Qtot,l,a.) 



1 /wh + Qtot 



+ Otot /?2-56fcpwKff/?2' 
Qtot /?fl'sinc(6'o) sin(a) 



fc2 (7/wh + Qtot i?2) 



(24) 
(25) 



One useful fact, made apparent by these equations, is 
that one can tune the parameter p independently of a 
by varying the wheel's inclination angle a. Furthermore, 
using these equation, we can determine the portion of 
parameter space that can be accessed with our imple- 
mentation of the water wheel. The accessible region's 
approximate upper boundary is shown in Fig. 4 by the 
black solid line. 

Although it is experimentally not possible to resolve 
the fine fractal structure seen in Fig. 4 because of ex- 
perimental limitations such as noise and parameter drift, 
we do find that our water wheel does qualitatively con- 
form to the theoretical predictions. For small p or cr, 
it shows steady state behavior of the angular velocity, 
corresponding to either stationary behavior or rotation 
with constant speed. For larger parameter values we find 
chaos or periodic oscillations. Furthermore, for fixed val- 
ues of (7 (cr ~ 3), we find that the wheel is in a chaotic 
state for small angles a (small p) and exhibits periodic 
behavior for larger angles a (large p). 

As an example of the obtainable data, we display in 
Fig. 5 two experimental time series, one periodic (Fig. 5a) 
and one chaotic (Fig. 5b). For each run several hours of 
data are taken, recording the wheel's angle as a function 
of time. The first roughly one hour is discarded as tran- 
sient and the rest is used for further time series analysis. 
Next, we low-pass filter and take the first derivative of 
the data in the Fourier-domain. Working in the Fourier 
domain is an efficient way of suppressing high frequency 
noise that would otherwise dominate the derivative. It is 
a legitimate way of proceeding because the data is highly 
oversampled (on the order of 100 points per oscillation 
period) and the Fourier spectrum has a dominant (but 
broadened) peak even for chaotic time series, which en- 
ables us to chose a filter that has essentially no effect at 
dynamically relevant time-scales. For example, the angle 
data corresponding to the chaotic time series in Fig. 5b 
is peaked around 0.07 Hz and the filter cutoff was chosen 
at ~ 0.6 Hz. Having computed the derivative, i.e. the 
angular velocity w, we obtain x by using the coordinate 
transformation (21), x = Uj/k. 

The upper (black) traces shown in Fig. 5 are ~7 min 
windows of data obtained after such processing. The time 
trace in Fig. 5a is clearly periodic. The small fluctuations 
of the amplitude are a result of the unavoidable imper- 
fections of the experiment. These become even more ap- 
parent in the two dimensional time-delay embedding^^'^^ 
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of the same data that is shown in Fig. 5c. For an ideal 
noise-free experiment exhibiting periodic oscillations, a 
single closed loop (a limit cycle attractor) should be seen, 
whereas our actual data does result in a smeared-out 
loop. Nevertheless, the correspondence between exper- 
iment and model is excellent, as can be seen in Fig. 5a 
by comparing the upper and the lower trace, which are, 
respectively, the experimental data and a numerical so- 
lution of the Lorenz equations [Eq. (22)] . 

The time scries in Fig. 5b is clearly non-periodic. A 
time-delay embedding of the same data is depicted in 
Fig. 5d and shows a complicated (potentially fractal) at- 
tractor. Together this indicates that the wheel oscillates 
chaotically. For comparison, we also display in Fig. 5b 
the result of a simulation of the Lorenz equations with 
the corresponding parameter values of p and cr, which lie 
in the chaotic regime of the model (see square in Fig. 4) . 
Again, we find a convincing qualitative correspondence 
between experiment and model in the sense that ampli- 
tude, dominant frequency, and general character of the 
solutions match. Since the solutions are chaotic they are, 
of course, not expected to be identical. 

Thus, qualitatively we find a good correspondence of 
the observed water wheel behavior and the predictions of 
the Lorenz equations. In order to obtain a more precise 
statement about the correspondence of model and exper- 
iment, we utilize another fascinating aspect of nonlinear 
dynamics: chaos synchronization. 



V. CHAOS SYNCHRONIZATION 

Though the synchronization of clocks (periodic sys- 
tems) has been studied with great care over centuries,^"* 
the discovery**^ that two chaotic oscillators could syn- 
chronize came as a surprise. Synchronization was unex- 
pected because sensitivity to small perturbations means 
that the solutions of two chaotic systems tend to diverge. 
Yet, it was found that, if two chaotic systems are cou- 
pled to each other in a suitable way, it is possible for one 
chaotic system to follow exactly the; time evolution of an 
identical second chaotic system, even when the chaotic 
systems start from very different initial conditions. In- 
deed, for the Lorenz system one can even prove that this 
is the case. 

The phenomenon of chaos synchronization suggests the 
following test of the validity of the Lorenz equations as a 
description of the water wheel. Presume that the water 
wheel dynamics is described by the Lorenz equations. 
Then a computer model has to synchronize to the water 
wheel for a suitably chosen coupling (e.g. for sufficient 
coupling strengths). That is, if we measure the angular 
velocity of the water wheel as a function of time, thereby 
determining x{t) in Eq. (22), and use this data as input 



to a numerical Lorenz model given by 

a;m = o- (ym -x^)-K {x^ - x) 

(26) 

where K sets a sufficiently large coupling strength, the 
model output Xai{t) has to converge to x{t), i.e. x^ — >■ x 
as t oo. If no convergence is observed, then the water 
wheel dynamics is not accurately modeled by the Lorenz 
equations. 

Before describing how our wheel performs under this 
test, let us discuss the proof of synchronization. 



A. Proving Synchronization 

First note that exactly synchronized solutions of the 
coupled system, Eq. (22) and Eq. (26), correspond to a 
zero error-vector, e{t) = 0, where the components of e 
are defined as (e^, ey,ez) = {x^ — x, — y,Zy^ — z). It is 
therefore convenient to consider the evolution of e, which 
is described by a non-autonomous ordinary differential 
equation of the form e = F{e,t), given by 

Cx — (Cy ^a:) ^ 

ey = pCx- ey - CxE:, - x{t) ez - z{t) (27) 
ez = exey-ez+ x{t) Cy + y{t) e^. 

Solutions e{t) can be thought of as trajectories in a three- 
dimensional phase space. To prove synchronization, we 
need to show that for any initial condition the corre- 
sponding solution trajectory asymptotically approaches 
the origin of this phase space, e — > as t — )■ oo. Since the 
solutions of the Lorenz equations (22) enter in Eq. (27), 
it is useful to establish bounds for the Lorenz solutions. 
Bounds can be obtained by constructing a trapping re- 
gion, that is, a closed region in the phase space of Eq. (22) 
with the property that solutions starting on the inside of 
this region never leave it and solutions starting on the 
outside enter it after some time. For the Lorenz equations 
with b = 1 [Eq. (22)] the following inequalities define a 
trapping region^^ 

x^ <4ap y"^ < p^ < z < 2p, (28) 

for positive real-valued parameters a and p. 

Next, we utilize a function that decreases along trajec- 
tories of system (27), a so-called Lyapunov function. As 
a Lyapunov function candidate consider the scalar func- 
tion 

p2 _(_ p2 _|_ p2 

F(e) = + (29) 

which is positive definite for e ^ 0, zero for e = 0, and 
monotonically increasing as a function of |e|. If we can 
show that F < for all e 7^ 0, then all trajectories flow 
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"downhill" toward the origin and e = is globally asymp- 
totically stable. Using Eq. (27) we obtain 



^(g) — Cx ~t~ ~t~ 



for the time derivative of V. Utilizing Eq. (28), a lower 

bound for the expression in the square brackets is found 
in terms of the time-independent parameters p and a, 



K + a- 



{a + p- zf y2 



- V >K + a- 
4 



{a + pf ^2 



4 ' 
(31) 



which, in turn, shows that for sufficiently large coupling 
strengths K, 



AK>a'^ + 2(7 p + 2p2 - 4cr, 



(32) 



the square-bracket term is positive definite. Then it fol- 
lows from Eq. (30) that 



F(e)<0 Ve^O. 



(33) 



This concludes the proof. 

It is important to emphasize that above proof and, in 
particular, condition (32) are sufficient conditions that 
guarantee that two Lorenz systems coupled unidirection- 
ally, as described by Eq. (26), will identically synchronize. 
Generally, the thus obtained value of -fC is a conservative 
estimate. Systems often synchronize for much smaller 
coupling strengths. It should also be noted that, in the 
K ^ oo limit, the model output will converge to any 
given input x, independent of the form of the dynamical 
system that generated x. In other words, the conver- 
gence of Xm to input X does not prove the correctness of 
the Lorenz model. However, if x^ does not converge to a 
(noise-free) input x for the sufficient coupling strength K 
derived above, x is definitely not generated by a system 
that evolves according to the Lorenz equations. 



B. Synchronizing the Model and the Wheel 

We now discuss results of coupling the chaotic data, 
part of which is shown in Fig. 5b, to the Lorenz model 
given by Eq. (26). When utilizing the provably sufficient 
coupling strength oi K = 2260 suggested by Eq. (32), 
we do find high-quality synchronization. What is more, 
we find that the model synchronizes to the data for much 
lower coupling strengths. An example is shown in Fig. 6a, 
where we display experimental data and model output for 
a coupling strength oi K = 100 and find that they are 
indistinguishable to the eye. The difference = — 
shown in Fig. 6b, has a magnitude that is about 1% of the 



a) 20 

10 

H 
-10 
-20 





20 30 

Time (arb.units) 

FIG. 6: (a) Experimental time-series (solid line) and model 
output (crosses), (b) Difference between model and experi- 
ment. 



oscillation amplitude (note the difference in scale) and is 
of much higher frequency than the natural timescale of 
the chaos, consistent with the interpretation that noise 
in the data is the main source of the remaining error. 

The fact that synchronization takes place for the suf- 
ficient coupling strength of the proof and, as it turns 
out, for miich smaller coupling strengths means that our 
wheel passes the test, it is not inconsistent with the 
Lorenz model. Indeed, within the limitations of the ex- 
periment, the presented data provides strong evidence 
that the Lorenz model is a good description of the water 
wheel dynamics. 



VI. DISCUSSION 

We have provided details on our inexpensive and easy 
to implement realization of the Malkus water wheel and 

demonstrated that it produces chaotic and periodic be- 
havior as predicted theoretically. We also showed that a 
numeric model implementing the Lorenz equations syn- 
chronizes to the wheel's chaotic motion with high fidelity. 
We discussed how chaos synchronization provides a test 
that, when not satisfied, can falsify the claim that the 
wheel is described by the Lorenz model. Within the lim- 
itations of an experimental system, our wheel passes this 
test. 

No model of a real system can be perfect and noise 
is present in any measurement. It is therefore not sur- 
prising that we find small differences between the data 
and model. The interesting question and open challenge 
is to develop tools that allow one to determine the main 
source of the remaining error. 

Perfect chaos synchronization between two (noise-free) 
chaotic systems is only achieved if the master system (the 
wheel) and the driven system (the simulation implement- 
ing Eq. (26)) are exactly identical. One might therefore 
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suspect that the parameters that were used in simula- 
tions were not exactly those corresponding to the wheel's 
operating condition because there undoubtably is uncer- 
tainty in the experimentally determined parameters that 
are used to map the measured angle 9{t) and wheel op- 
erating conditions onto the variables and parameters of 
the Lorenz system. Stated differently, one is confronted 
with the problem of finding optimal parameter values for 
a chaotic system based solely on ones knowledge of a 
single scalar time series, in this case 9{t). This is a diffi- 
cult problem, in general, due to the sensitive dependence 
of chaotic systems on the parameters and initial condi- 
tions of the unknown state variables {y and z) . Recently 
proposed techniques, that are inspired by control engi- 
neering approaches, are in some cases able to solve this 
problem. ^^"^^ Utilizing these ideas, we designed an adap- 
tive observer^^ that generates for a given time series an 
(optimal) estimate of the corresponding model parame- 
ter values. It is these optimized values that were used 
for all simulation results shown in this paper. A detailed 
discussion of the adaptive observer technique is beyond 
the scope of this paper. However, we note that, quite 
reassuringly, the optimized parameter values are in close 
correspondence to the estimates based on experimental 
measurements. For the chaotic time series the optimized 
values are a « 2.5 and p « 66 as compared to the exper- 
imentally determined values of cr w 2.7 and p « 69. As a 
matter of fact, given the measurement uncertainties we 
find this close correspondence rather surprising. 

In summary, of the three potential causes for the small 
differences between the data and model, namely noise, 
parameter mismatches, and structural insufficiencies of 
the model, we can exclude parameter mismatches. We 
know for certain that the data is noisy and therefore con- 
tributes to the observed differences. We cannot with cer- 
tainty exclude model insufficiencies, but if they exist their 
effect is small. Thus, we think that the wheel described 
in this paper comes close to an ideal one, one whose dy- 
namics is exactly described by the Lorenz equations. It 
can therefore be used to explore other interesting aspects 
of the Lorenz system, such as its bifurcations. 
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Appendix A: Numeric calculation of the Lyapunov 
Exponents 

For completeness we will recall in this section the def- 
inition of Lyapunov exponents and provide some details 
on how to compute them numerically. Some reviews on 
this topic are found in Wolf et al.,^°, Geist et al.,^^ and 
Souza-Machado et al.^^ 

The object of study are nonlinear ordinary differential 
equations, 

X = f(x), x(0) = xo, (Al) 

where x(t) = (a;i(t), .X2(t), . . . , x„(t)) G M", xq is the 
initial condition, and f is an n-dimensional vector field. 
To define Lyapunov exponents associated with solutions 
x(t), one needs to take into account that the rate of 
separation of two trajectories in phase space can be dif- 
ferent for different orientations of the initial separation 
vector. Therefore, one has to consider a spectrum of n 
Lyapunov exponents for an n dimensional phase space, 
{Al, A2, . . . , A„}. Since the Lyapunov spectrum charac- 
terizes the evolution of infinitesimal perturbations, one 
can utilize a linearization around a fiduciary trajectory 
via solutions to the matrix differential equation 

Y = J Y; Y(xo, t = 0) = l. (A2) 

In this equation, J is the time-dependent Jacobian ma- 
trix with elements Jij[x(t)] = (9/i/9.Tj)|x=x(/,) that are 
the partial derivatives of the vector field f in Eq. (Al) 
evaluated along the fiduciary trajectory x(f). With the 
identity matrix as initial condition the n x n matrix Y 
gives then the complete linearized flow map with re- 
spect to the standard orthonormal basis {ei, . . . , e„}. In 
other words, Y(xo,t) describes the evolution in time of 
both the magnitude and the orientation in phase space 
of any initial infinitesimal perturbation from xq because 
Sx{t) = Y(xo,t)(5xo. The Lyapunov exponents Aj are 
defined by the logarithms of the (real) eigenvalues fii of 
the positive and symmetric matrix 

A = hm [yT(xo, t) Y(xo, t)] , (A3) 

where Y'^{xo,t) denotes the transpose of Y(xo,f). This 
implies that for every initial condition xq there exists an 
orthonormal set of initial perturbation vectors 5vi such 
that 

Ai = lim ^ ln|Y(xo,t)5vi| , i = l,2,...,n. (A4) 

It was shown by Oseledets'^^ that the limits on the right- 
hand side of Eq. (A3) and Eq. (A4) exist for almost every 
initial condition xq and it has been argued^^ that for er- 
godic systems the values of the Lyapunov exponents {A^} 
do not depend on the initial conditions (up to a measure 
zero in phase space). Thus, the Lyapunov exponents are 
global properties of the attractor of the dynamical sys- 
tem. 
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The definition of the Lyapunov exponents might sug- 
gest to integrate the iinearized equation (A2) along with 
the nonhnear differential equation (Al). This is not fea- 
sible because the exponential growth and decay rates of 
initial perturbation vectors quickly exceed the abilities 
of numerical number representation in a computer. Fur- 
thermore, perturbation vectors quickly align with the di- 
rection of maximal growth making it impossible to de- 
termine any but the largest Lyapunov exponent. To 
compute the Lyapunov spectrum, the most commonly 
used algorithms are based on a step- wise procedure where 
Eq. (A2) is integrated for short intervals of time and the 
vectors spanning the phase space volume are reorthonor- 
malized at each step. 

To be precise, we compute the largest k Lyapunov ex- 
ponents by considering the evolution of the volume 14 
of fc-dimensional parallelepipeds in n-dimensional phase 
space. It can be shown that 14 is given in terms of the 
sum of the largest k Lyapunov exponents^^'^^ via 

k 

V Afc = lim I In [14] (1 < fc < n) (A5) 

i=l 

when the Aj form a monotonically decreasing sequence. 
Consider therefore an initial fc-dimensional hypercube 
centered on the initial point Xq, where the cube is de- 
fined via an orthogonal n x k matrix Pq = (6i, . . . , 6k), 
the columns of which are k (randomly chosen) orthonor- 
mal n-dimensional vectors 6^ forming the cube's axes. 
Under the action of the flow the cube will deform into an 
parallelepiped P(t) = Y(xo, i)P(O) with volume I4. The 
volumes 14 can be computed via the QR decomposition 
that factorizes a n x k matrix P 

P = QR (A6) 

into the product of an orthogonal nx k matrix Q and an 
upper triangular k x k matrix R with positive diagonal 
elements Ru. The volume of P is 

k 

14 = I det P| = I det R| = JJ Ru (A7) 

i=l 

Substituting this expression into Eq. (A5) for I4 for all 
k with 1 < A; < n, we find that 

Xi = lim I In [Rii{t)] {i = l,...,k) (A8) 

We can calculate the time average of the Ru in 
Eq. (A8) in a stable manner. To do so, first note that 
the evolution in time of P is given by 

f = ^P(0) = JPW. (AC,) 

The advantage of considering the time evolution of P is 

that the determination of the largest k Lypunov expo- 
nents involves the integration of just n{l + k) differential 



equations instead of the n (1+n) equations that would be 
necessary if one considered the evolution of Y. The com- 
putation proceeds stepwise, where one integrates over a 
sufficiently small time-interval, Atj = tj — tj-i, the cou- 
pled differential equations 

H-x: 

-=f(x), X(t+_,)=x(i7_j 

jp tj-i <t<tj. 

— =J[x(t)]P, P(t+_,) = Q,_i, 

(AlO) 

The initial conditions for the first step are x(0) = xg 
and P(0) = Pq. The size of the time step Atj has to be 
chosen such that P{t) remains well conditioned. Time 
steps on the order of a typical oscillation period are often 
suggested. Next, the QR decomposition of the matrix 

= Pit 3), 

P, = Q,-R,-, (All) 

is performed to reorthonormalize the phase space vol- 
ume. The orthogonal matrix Qj is used to initialize the 
integration of the next step and the diagonal matrix ele- 
ments of Rj are accumulated to compute the Lyapunov 
exponents. Since the fiow matrix P = P^ • • • Po can be 
expressed as the product of the matrices P^ computed at 
successive points along the orbit, x{tj), above procedure 
implies that 

P = Qm Rm Rm-l • • • Rl- (-^12) 



Thus, the Lyapunov exponents are 




FIG. 7: Geometric illustration of the QR decomposition 
scheme for computing Lyapunov exponents. 

We implement the above procedure using an adaptive 
stepsize integrator (the DDRIV fortran-integrator) and 
the QR decomposition routine contained in the LAPACK 
library. 

In general the convergence to the limit in Eq. (A13) 
is very slow and it is useful to have checks on how well 
the exponents have been approximated. To do so, we 
note that for a dissipative systems, such as the chaotic 
Lorenz system, the total phase space volume contracts 
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(exponentially in time on average), implying that the sum 
of all Lyapunov exponents, 

3 1 1 

y A, = lim - InldetYl = lim - \n[V3(t)] , (A14) 



i=l 



must be negative. For the Lorcnz system with & = 1 
[Eq. (22)] with corresponding linearized flow equation 



Y = JY 




(A15) 



one can calculate the phase space volume contraction rate 
explicitly by applying to equation Eq. (A15) Liouville's 
formula 

det Y(xo,i) = det Y(xo,0) exp^^ tr[J(0] d^^. 

(A16) 



Since dctY(xo,0) = det 1 
cobian is a constant, tr[J] 
that 



1 and the trace of the Ja- 
-a - 2, Eq. (A14) implies 



J2^i= lim ] f tr [J(0] dC = -<7 - 2. (A17) 

^ t^cx) t Jo 



This provides a cross-check for convergence. A second 
check is obtained by noting that perturbations along a 
trajectory will neither grow nor shrink exponentially in 
time, implying that there is a zero Lyapunov exponent. 
This zero exponent exists whenever the attractor is not 
a fixed point, e.g. for periodic or chaotic dynamics. 

We utilize these checks in our numerical procedure by 
computing all three Lyapunov exponents for the Lorenz 
system and accepting the computed values whenever (a) 
the zero exponent, if it exists, has converged such that its 
magnitude is smaller than 10"^ and (b) the sum condi- 
tion Eq. (A17) is satisfied, -{2 + a)- Aj < IQ-^. 
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